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Abstract. An abstract should be given 



Exact analytic solutions are derived for radiation 
in time-dependent relativistic flows. The flows are 
spherically-symmetric homologous explosions or implo- 
sions of matter with a grey extinction coefficient. The so- 
lutions are suitable for testing numerical transfer codes, 
and this is illustrated for a fully relativistic Monte Carlo 
code. 
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1. Introduction 



methods:analytical 



In a recent paper (Lucy 2004), a Monte Carlo (MC) treat- 
ment accurate to 0(v/c) of the time-dependent transport 
of radiation in 3-D SNe is described and tested. A major 
concern in that paper is establishing the accuracy of the 
MC code. To this end, the 3-D code was applied to a 1-D 
problem that could be solved independently with conven- 
tional numerical methods. Specifically, the test problem 
was to compute the bolometric light curve of a spherical 
SN in which the transfer of UVOIR radiation is treated 
with a grey extinction coefficient. An independent ap- 
proach to this problem is provided by Castor's (1972) 
co-moving frame (cmf ) moment equations for spherically- 
symmetric flows. The resulting pair of partial differential 
equations (PDEs) were solved with the Henyey method. 

The final outcome of this test was entirely satisfac- 
tory: the mean difference between the two light curves 
is ^ O.Olmag. for elapsed times t from 10 to 50 days. 
Nevertheless, the differences were initially significant, and 
it was not clear which code was in error. This question 
was eventually answered by testing each code separately 
against an exact similarity solution of Castor's equations. 
The cause of the differences could then be traced to the 
poor spatial resolution of the MC code. 
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In the interest of concise presentation, the essential 
part played by this similarity solution was not described 
in the earlier paper (Lucy 2004). But subsequently, this 
solution was found to generalize to all orders of v/c. 
Accordingly, since time-dependent relativistic flows and 
the associated transfer problems are of interest for such 
phenomena as gamma-ray bursts and micro-quasars, this 
paper derives this solution (and variants thereof) and il- 
lustrates its use in testing a relativistic transfer code. 

2. Basic equations 

Similarity solutions will be sought for the frequency- 
integrated radiation field in a homologously expanding 
or contracting flow. The configuration is spherically- 
symmetric, and the matter has grey extinction (\) and 
integrated emissivity (77) coefficients that are isotropic in 
the cmf. 

2.1. The transfer equation 

Mihalas (1980) has derived the general time-dependent 
cmf transfer equation for spherically-symmetric flows with 
relativistic velocities. From this, he derives the zeroth and 
first frequency-integrated moment equations, whose ear- 
lier derivation by Prokof 'ev (1962) is acknowledged. When 
terms of 0(v 2 /c 2 ) and higher are neglected, Castor's 
(1972) moment equations are recovered. 

Following Mihalas (1980), the interaction coefficients 
(x, T]) and the radiation field (/, J, -ff, K) are expressed in 
the cmf but spacetime coordinates (r, t) and flow velocities 
(y) are measured in the rest frame (rf). Since radiation 
quantities always refer to the cmf, primes or suffixes to 
indicate this frame are omitted (cf. Mihalas 1980, Sect. 
III). 

In a homologous spherical flow, we have v — r/t. For 
t > 0, the flow is an explosion (v > 0) starting with infinite 
density at t = 0. For t < 0, the flow is an implosion (v < 0) 
leading to infinite density at t = 0. 
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When Eq. (2.12) of Mihalas (1980) is applied to the 
homologous flow of matter with a grey extinction coeffi- 
cient, the resulting frequency-integrated transfer equation 
is 

7/., n \dl u dl 7 ,„ n.dl 47 



= r]-xl 



(I) 



Here the dependent variable is the integrated specific in- 
tensity I{ji;(3,t), the radial coordinate r has been re- 
placed by (3 = v /c, the time derivative is Lagrangian, and 

7 = l/V / l^- 

2.2. Moment equations 

Moment equations can be derived from Eq. (I) as usual 
by multiplying by \[i m with m = 0, l, 2, ... and then inte- 
grating over direction cosine \x. Alternatively, Eqs. (2. 16) 
and (2.17) in Mihalas (1980) can be simplified to the case 
of homologous flow, as above for his transfer equation. 
The resulting equations are 

7_riJ jP_dH_ J_dH_ 27 

cdt + c dt + yet 8(3 + Pct [U+2 ' jJ) 



for the zeroth moment, and 

7 riff 7/3 dK J_dK_ 7 

c dt c dt yet d[3 pet 



= v-xJ 



{3K - J) 



(2) 



-XH 



(3) 



for the first moment. 

Following Mihalas (1980), Eddington's flux variable H 
is here preferred to the standard flux F = 4H. The inte- 
grated physical flux in the cmf is therefore T = 4wH . 

Comparing Eqs. (2) and (3) with Eqs. (2.16) and (2.17) 
of Mihalas (1980), we see that specializing to homologous 
flow has resulted in a huge simplification in the coefficients 
of the moments J, H and K . Most notably, the coefficient 
of K in the zeroth moment equation reduces to zero, so 
now only J and H appear in this equation. 

2.3. Separation of variables 

In order to construct simple solutions of these equations 
for use in testing computer codes, we first effect a separa- 
tion of variables. This is achieved by writing the specific 
intensity 



I(Mp,t)=h(»;p) {h/tf 



(4) 



where t\ is an arbitrary reference time, and the exponent 
p is unspecified. With this assumption, each term on the 
left-hand side of Eq. (1) scales as (ti/t) p+1 . Accordingly, 
we must assume that 



V (p,t)= m (p) (h/ty +1 

and 

x(P,t) = xi(P) {hit) 



(5) 



Not surprisingly, this necessary scaling of the extinction 
coefficient per unit volume implies that the ratio of 1/%, 
the local mean free path of a photon, to the radius of the 
configuration is time-independent. 

When Eqs. (4)- (6) are substituted in Eq. (1), the scal- 
ing factor cancels by construction, and we thus obtain the 
transfer equation 



M 9h ' 7 -(l-M a )^ + £[4-Ml+ </<>]/, 



+ 



7cii dp Pch y ^ ' d[i 

= Vi-Xi h (7) 

that determines the scale- free radiation field 7i(/z;/3). 

The moments J , H and K evidently also scale as 
(ti/t) p . The equations satisfied by the corresponding 
scale- free moments Ji(/3), and K\{P) are 

jcti dp Pct\ ct\ 

= Vi~ X1J1 



(8) 



for the zeroth moment, and 
1 dK 1 



jP R 
-fcti dp ct\ pct\ 



(3Jfi- Ji) + ^(4-p)ffj 

= -X1H1 (9) 

for the first moment. 

Eqs. (7)-(9) are the basic equations of this investiga- 
tion. 

3. Testing a Monte Carlo code 

If MC techniques are used directly to simulate the physics 
of radiation transport, then the MC quanta are photons 
and convergence to the solution of the Radiative Transfer 
Equation (RTE) requires that Af — > 00, where Af is the 
number of photons whose interaction histories are fol- 
lowed. In such a code, in addition to crossings of bound- 
aries, MC quanta are created spontaneously within the 
computational domain D by sampling the thermal emis- 
sivity and may subsequently be destroyed within D by 
absorption. Because the RTE is directly simulated, test- 
ing such MC codes is not fundamentally different from 
testing a conventional numerical solution of the RTE. 

However, for tranfer problems involving interactions 
between radiation and the internal energy states of mat- 
ter, there are advantages in taking the MC quanta to be 
indestructible and indivisible energy {£-) packets (Lucy 
2003, and references therein). In such a code, in addition 
to crossings of boundaries, ^-packets are created spon- 
taneously within D by sampling the net emissivity (i.e., 
emission minus absorption) but then, though the nature 
of the contained energy may change, they are not subse- 
quently destroyed within D by absorption. 

3.1. Moment solution 

Let the integrated net emissivity per unit volume at time 
t\ be 47r?7i, where 



(6) fji{p) = Tii - X1J1 



(10) 
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This quantity creates radiant energy within the configu- 
ration by a physical mechanism that need not be specified 
for test problems. 

If we replace the right-hand side of Eq. (8) by rji and 
eliminate J\ from the left-hand side by setting p = 4, 
the result is an ordinary differential equation (ODE) for 
Hi(J3), 

( ^l + 2 -f-{l-2(i 2 )H l ^ 1 ct 1 f ll (11) 

This equation can be solved analytically with the inte- 
grating factor (3 2 /-f 2 . The solution satisfying the boundary 
condition i?i(0) = is 

#i(/3) = ch 4 / m(b) &V1-6 2 db (12) 
P Jo 

Note that this formula is exact. In transfer theory, ana- 
lytic formulae for moments are typically not exact because 
they are derived with Eddington's closure approximation 
K = J/3. Here this is not necessary: K and J drop out of 
the zeroth moment equation because of the assumptions 
of homologous flow and scaling exponent p = 4. 

A second point to note is that -ffi(/3) has been de- 
rived without specifying the scale-free extinction coeffi- 
cient. Accordingly, Eq. (12) is valid for arbitrary Xi(P)- 

It is of interest to note that the scaling p — 4 arises 
naturally in the limiting case of completely opaque matter 
within which there is no net emissivity. In this case, the 
right-hand side of Eq. (2) is zero, and the cmf flux H = 
since radiation is position-coupled to matter. The solution 
of Eq. (2) is then such that J oc t~ 4 , corresponding to 
adiabatic evolution of the radiation energy density present 
initially. In contrast, in the solution derived here, this same 
scaling is maintained because every layer's losses due to 
flux divergence is exactly replaced by the net emissivity in 
that layer. Moreover, the similarity solution represents the 
state reached when initial conditions have been erased. 

3.2. A particular case 

The simplest case for testing a MC code is when ryi(/3) is 
independent of j3. This also simplifies the evaluation of the 
exact Hi((3) since the integral in Eq. (12) is then analytic. 
The result is 

H 1 {(3) = \ct 1 f,i Xr [^n- 1 /3--(l-2/3 2 )] (13) 

Thus we obtain an exact closed-form expression for the 
cmf flux in a particular time-dependent relativistic flow. 

The corresponding moment J\ (j3) cannot be obtained 
exactly. But an approximate formula can be derived from 
Eqs. (9) and (13) with the help of Eddington's closure ap- 
proximation and surface boundary condition. The result- 
ing formula for Ji can then be substituted in Eq. (10) to 
obtain an approximate formula for the conventional emis- 
sivity rji. Details are omitted. 



3.3. Monte Carlo calculation 

In order to illustrate how similarity solutions can be 
used to test codes, we now briefly report MC calcula- 
tions for relativistic homologous flows. The MC code is 
a spherically-symmetric and fully relativistic version of 
the 3-D code described recently (Lucy 2004). As in that 
code, the MC quanta are indestructible ^-packets. The 
calculations start at reference time t\ with no ^-packets 
present. But as time advances, ^-packets spontaneously 
appear in accordance with the net emissivity fji(ti/t) 5 and 
then propagate through the configuration interacting with 
matter in accordance with extinction coefficient \i(ti/t). 

We choose to create equal numbers of ^-packets in 
equal intervals of logt and having cmf energies e(t) that 
are independent of (3. The energy dE created in the cmf 
within the space-time element dV dt is dE = 4irf)dV dt. 
But since dV dt is invariant under the Lorentz trans- 
formation, we may regard dV dt as referring to the rf. 
Accordingly, for the particular case (Sect. 3.2) where fj 
docs not depend on /3, the total cmf energy created in the 
rf interval dt is AwrjV dt, where V(t) = Vit^^/hf is the 
rf volume at time t. Thus, if dAf/dint gives the packet 
creation rate, and we set p — 4, then the packets' cmf 
energies are 

^ = ^■(7)®" <14) 

With f] independent of /?, dE oc dV at fixed rf t. 
Accordingly, since we choose e to be independent of /?, 
equal numbers of packets are created in equal rf volume 
elements. The appropriate MC sampling algorithm for the 
initial rf radius of a packet created at t is r = i?* \fz, where 
i?* = c(3*t and z is a random number from (0, 1). Here /3* 
is the value of v/c at the surface. 

The transport of these packets through the configura- 
tion is as described in the earlier paper (Lucy 2004) except 
that now the factor 7 in the Doppler formula is restored 
to make the code fully relativistic. The cmf energies of the 
packets that escape from the surface in each rf time step 
At are summed and then divided by At/7* to obtain an 
estimate of the cmf luminosity L(t) 

Monte Carlo simulations for two strongly relativistic 
explosions are plotted in Fig. 1. For these simulations, the 
rf time steps are Alogt = 0.02, during each of which Af = 
1000 additional ^-packets are created with cmf energies 
given by Eq. (14). 

The evolution of the cmf luminosity is shown start- 
ing at t\, together with the predicted result L = A-kR 2 x 
4ttFi(/3») x (ti/t) 4 , with #i(/3*) from Eq. (13). One so- 
lution is the free streaming limit (xi = 0) and the second 
(xi — 3) has photon mean free paths = i?*(i)/3. In both 
cases, the MC solutions tend asymptotically to the sim- 
ilarity solution, the convergence being somewhat slower 
for xi = 3 because of the longer residence times of the 
^-packets. Since agreement is achieved in both cases, the 
prediction, for p = 4, that the cmf flux H is independent 
of xi is confirmed. 
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Fig. 1. Explosions. Comparison of MC calulations (cir- 
cles) of the cmf luminosity L(t) with predictions of similar- 
ity theory (straight lines). The surface velocity v* = 0.9c. 
The unit for the indicated scale-free extinction coefficients 
Xi is l/Ri, where Ri — v*ti. The unit of luminosity is 
4/37ri?f x 471771. The luminosites for xi = 3 have been 
increased by 1.5dex. 

3.4. Relativistic implosions 

Although relatively trivial, it is of interest to illustrate 
the application to implosions. Accordingly, Fig. 2 shows 
the same two test problems but with the signs of t\ and 
/3* reversed. Again, after a transition period during which 
the internal radiation field is established, the MC solutions 
converge to the similarity solution. 

This possibility of testing numerical treatments of ra- 
diative transfer in relativistic implosions is possibly rel- 
evant for codes that simulate laser-driven implosions in 
support of the quest for fusion by inertial confinement. 

Accurate numerical solutions for relativistic inflows 
have also been computed by Yin & Miller (1995), who 
stress the dramatic effects that can arise from photon trap- 
ping. But their solutions are only for stationary flows. 

4. Testing a transfer code 

Apart from the degenerate case Xi = 0, the exact solution 
of Sect. 3 is not appropriate for the precision testing of 
a conventional relativistic transfer code since the implied 
scale-free emissivity r\\ could only be determined approxi- 
mately. Accordingly, we now seek exact solutions when r\\ 
is specified rather than f\\. 



4.1. Solution along characteristics 

As did Mihalas (1980, Sect. Illb) for his general equa- 
tion, we construct characteristics for Eq.(7) such that the 
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Fig. 2. Implosions. Same as Fig. 1 but with signs reversed 
for [3* and t\. The luminosites for xi — 3 have been in- 
creased by 1.5dex. 

partial differential operator becomes a perfect differential. 
Thus Eq. (7) becomes 



1 dh 
+ 

cti ds 



Xi 



(p-4)7 + p7/3/z 



cti 



h = m 



(15) 



where ctis is distance along a charactcristc defined by the 
equations 



<W_^!i and jv_ = 2n 

ds 7 ds (3 



l' 



(16) 



Dividing these equations, we find that ^f3^J(l — y?) is 
constant along characteristics. Accordingly, the family of 
characteristics is given by 



7 /V(l-/i 2 )=7*/?*V / (l-^)=a 



(17) 



where — 1 < /i* < is a characteristic's direction cosine at 
its entry point into the configuration. Integrations of Eq. 
(15) proceed inwardly along such characteritics starting 
at s — with the boundary condition h(fi*; /?*) = and 
continuing until the characteristic emerges at the surface 
with \i= 

In Fig. 3, the characteristics are plotted for various 
values of /j» when /?* = 0.9. Note that, to obtain a charac- 
teristic that penetrates close to the centre, the parameter 
/i» must closely approach —1. The point of closest ap- 
proach (ji = 0) is at 7 = Vl + a 2 = 70 , or, equivalently, 
at (3 = a/70 = (3o- 

Parenthetically, we note that these analytic character- 
istic curves provide a further powerful test of the relativis- 
tic MC code described in Sect. 3.3. In the free streaming 
case, each ^-packet follows its appropriate characteristic 
to high precision, as it should. At creation, an ^-packet's 
initial (3 and fx determine its invariant 7/3 -^/(l — fi 2 ) and 
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Fig. 3. Characteristic trajectories /u(/3;/i*) for Eq. (7) 
when /?* = 0.9. The curves are labelled with the value 
of the direction cosine at the surface for the inward 
ray. 

hence /i» from Eq. (17). The packet then propagates along 
this characteristic in the direction of increasing s until it 
escapes at the surface (3 — /3* with /i = 

The simplicity of the characteristics for homologous 
flow allows the dimensionless arc length s to be evaluated 
analytically as a function of (3. Integrating the first mem- 
ber of Eq. (16) after eliminating /i with Eq. (17), we find 
that, along the inwardly-directed (^ < 0) segment of a 
characteristic, 



«(/3;m*) = — 

7o 



sin ( — ) - sin ( — ) 

7 7* 



(18) 



The corresponding formula after the point of closest ap- 
proach (u > 0) is 



s(/?;/x*) = 



1 

7o 



n — sin ( — ) — sin ( — ) 

7 7* 



(19) 



Because the function s((3; /z*) is readily inverted to give 
/3(s;/«*), the remaining quantities 7, /x, 771 (/?) and Xi(/3) 
can likewise be transformed into functions of s along the 
characteristic defined by /z* for the given /3*; and this re- 
mark therefore applies also to the coefficient of Ji in Eq. 
(15). 

Let us now define the effective extinction coefficient 
along a characteristic to be 



xi(s;a**) = xi - 



(p-4)>y + pjf3u 
cti 



with corresponding effective optical depth 

fi(s;^») = cti / xids 
Jo 



(20) 



(21) 



In terms of these quantities, the formal solution of Eq. 
(15) subject to the boundary condition / = at s = is 



Ji(s;/i.) =cti / Vi(s')e f(s 



(22) 



With straightforward numerical integrations, Eqs. (18) - 
(22) allow the intensity I(s;u*) to be determined as a 
function s along the characteristic defined by /i» . By vary- 
ing u* from —1 to 0, we can thus determine / throughout 
the (/x, /3)-plane. 

This reduction of the problem to solution with the for- 
mal integral strongly suggests that these results could also 
be obtained with the method developed by Baschek et al. 
(1997). Although their paper is restricted to stationary 
relativistic flows, they note that time-dependent problems 
can be treated. 

4.2. Equivalent static medium 

The above analysis shows that, with the scaling assump- 
tions for n and % given in Eqs. (5) and (6), the time- 
dependent relativistic transfer problem for homologous 
spherical flow, reduces asymptotically to a tranfer prob- 
lem in a static medium. For explosions, the solution tends 
to this asymptote as t — > 00; for implosions, the limit 
is t — > 0. For free-streaming, the similarity solution is 
achieved after all light signals emitted from within the 
configuration at t = t\ have escaped. 

Because light travel-time and relativistic effects are 
absent in a static medium, these effects reappear in Eq. 
(20) as corrections to the extinction coefficient. The term 
oc (p — 4), which survives in the limit j3 — * 0, represents 
the combination of finite propagation speed and the time- 
dependence of the emissivity. The second correction term, 
which — > as j3 — > 0, represents relativistic effects. 

Note that these corrections can give xi < 0, which 
implies that f 1 is not necessarily a monotonically increas- 
ing function of s. Accordingly, fi is not appropriate as an 
alternative independent variable for Eq. (15). 

4.3. A particular case: moment solution 

Although the analysis of Sect. 4.1 provides a complete 
solution for the scale-free radiation field without assump- 
tions about n\{j3) or Xi(/3)> the answer is not in closed 
analytic form as was the earlier result (Eq. [13]) for the 
cmf flux H\ when the net emissivity fj is independent of fi. 
Interestingly, an analogous result can be constructed when 
the conventional emissivity n is the quantity specified. 
The coefficient of Ji in Eq. (8) is zero if we set 

7 



Xi(/3) = (p-4) 



ct-i 



(23) 



Thus p > 4 is necessary for xi > 0, while p = 4 gives 
free-streaming radiation. 

With xi given by Eq. (23), Eq. (8) simplifies to 



~dfT + l3 



(2-p/3 2 )H 1 = cti 7 Vi 



(24) 



6 



L.B.Lucy: Similarity solutions 



This equation has integrating factor 1 ' /j p 2 , and so the 
solution satisfying the boundary condition Hi (0) = is 



7 



p-2 rP 



a 2 



m (b) b\\-b 2 y 



db 



(25) 



A simple case for which this integral is analytic is obtained 
by assuming that rji is independent of (3 and that p = 5. 
The result is 



1 

15 



ch m 7 3 /3 (5 - 3/3 2 ) 



(26) 



With this exact formula, a code that solves relativistic 
transfer problems with conventional techniques can be 
subjected to tests similar to those described in Sect. 3.3 
for a MC code based on ^-packets. 

4.4. A particular case: complete solution 

The particular case of Sect. 4.2 is remarkable in that the 
complete solution can be obtained in closed form. The 
starting point is Eq. (15) with xi from Eq. (23). The in- 
dependent variable s is conveniently transformed to (3 us- 
ing the first member of Eq. (16). The resulting transfer 
equation is 

^-p 1 *(3I l =ct m ^ (27) 
dp [i 

where, from Eq. (17), 

m= ±J\/P 2 ~P 2 o (28) 

Here the negative root applies along the inward segment 
of a characteristic as [3 decreases from /?* to [3$ - see Fig. 
3. Thereafter, the positive root applies. 

The integrating factor for Eq. (27) is j~ p . Accordingly, 
substituting for /i, imposing boundary condition I~ (/?*) = 
0, setting p = 5, and assuming that r/i is independent of (3, 
we find that the intensity of inwardly-directed radiation is 



r(f3;n*)=ct im — I \ . ; „ bdb 



70 J p 

which simplifies to 

7 5 
7o 



where 
with 



(29) 



(30) 



(31) 



(32) 



The outwardly-directed intensity is obtained similarly. 
In this case, the boundary condition is / + (/?o;/i*) = 
J~(/?o;M*), and the solution is 



7 
7o 



(33) 



Emergent Intensities 

Homologous 
explosions 




Fig. 4. Emergent intensities in the co-moving frame for 
the indicated values of /3*. The exponent p = 5, and xi 1S 
from Eq. (24). The formula for Xi(P) m units of l/Ri is 
exhibited (cf. Figs.l and 2). The unit of intensity is ct\T]\. 



Eqs. (30) and (33) determine the intensity along the 
characteristic defined by //* . But these formulae can read- 
ily be used to calclulate I\{p,;(3), the intensity as a func- 
tion of jU at fixed (3. Given n and /?, Eq. (17) determines 
a, which in turn determines (3q, so that /3 = (3) and 
similarly of course for 70. Thus the right-hand sides of 
Eqs. (30) and (33) are now functions of \x and /3, and the 
characteristics are no longer relevant. 

This exact closed form analytic solution has been 
checked by substitution back into Eq. (7) using numer- 
ical differentation to evaluate the two partial derivatives. 

In Fig. 3, the emergent intensity given by Eq. (33) is 
plotted for various values of /?». This shows extremely 
strong forward peaking when /3* = 0.9, with huge inten- 
sities when [i <; 0.9. This is basically a light travel-time 
effect: with /?„ = 0.9 radiation emerging with \x <~ 1 in- 
cludes photons emitted shortly after the explosion when, 
since j oc 1/i 6 , the cmissivity was much higher than that 
now in the surface layers. 

One consequence of this forward peaking is poor ac- 
curacy for Eddington's approximations. Thus, at the sur- 
face of the 0* = 0.9 solution, K x j J\ = 0.821 and H x j J\ = 
0.895, as against Eddington's values of 1/3 and 1/2, repec- 
tively. 

4.5. A particular case: thermal emission 

Thus far the physical mechanisms responsible for extinc- 
tion and emission have not been specified. But let us now 
suppose that 



X = k + cr 



(34) 
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where k and a are grey absorption and scattering coef- 
ficients, respectively. The corresponding integrated emis- 
sivity coefficient is then 

r, = kB + nJ (35) 

where B is the integrated Planck function, and the scat- 
tering is isotropic in the cmf. 

Eqs. (34) and (35) are consistent with the scaling of 
Eqs. (5) and (6) if k and a are oc 1/t and if B oc l/t p . 
Given these scalings, the exact solution of Sect. 4.2 can 
now be applied as follows: we set p = 5 and choose the 
scale-free functions k\(fi) and 0i(/3) such that their sum 
satisfies Eq. (23). Then, since Ji(/3) can be computed from 
Eqs. (30) and (33), the implied scale- free function Bx{(3) 
can be derived from Eq. (35). Because the thermal emis- 
sivity is now known, this analytic solution can be used 
to test a code that incorporates an iteration procedure to 
detemine the scattering contribution to the emissivity. 

The net emissivity fji = r/i— XiJi = k(Bi — J\) is also 
determined by the above steps. Accordingly, this analytic 
solution can be used to test MC codes based on £ -packets. 
Moreover, this test is in principle more powerful than that 
provided by the flux moment solution of Sect. 3.2 since the 
angular distribution of the MC radiation field can now also 
be checked. 



higher order difference formula can be employed that uses 
the solutions at the two previous time steps (Lucy 2004). 
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5. Conclusion 

The aim of this paper has been to derive exact ana- 
lytic solutions in order to test radiative transfer codes 
for relativistic flows. Such solutions exist for spherically- 
symmetric homologous flows with power-law time depen- 
dencies for the grey extinction coefficient and the inte- 
grated emissivity. The exact solution for the integrated 
intensity derived in Sect. 4.4, being a function of three 
independent variables t,fi,(i and one parameter /?*, pro- 
vides an extraordinarily demanding and informative test 
for such codes. Moreover, its somewhat contrived deriva- 
tion based on a particular spatial variation of \ m no wav 
lessens its usefulness. 

In addition to deriving closed form analytic formulae 
for the intensity and flux in the cmf, an exact formula 
has been derived for the characteristics. This allows the 
kinematic and geometric aspects of relativistic transport 
codes to be tested independently of the treatments of ab- 
sorption and emission. For MC codes, this test is carried 
out by setting the extinction coefficient to zero and then 
checking that photon packets propagate along the known 
characteristics. 

In solving time-dependent transport problems with the 
RTE, the time derivatives are commonly approximated 
with a backward difference formula, thus using the so- 
lution at the previous time step. This introduces errors 
O(At) that might well accumulate as the integration pro- 
ceeds. With the availability of an exact time-dependent 
solution, the magnitude of the accumulated error can be 
determined. If the error accumulation is unacceptable, a 



